Numerical simulation study of the dynamical 
behavior of the Niedermayer algorithm 



D. Girardi, N. S. Branco 

E-mail: nsbranco@fisica.ufsc.br 

Departamento de Fi'sica, Universidade Federal de Santa Catarina, 88040-900, 
Florianopolis, SC, Brazil 

Abstract. 

We calculate the dynamic critical exponent for the Niedermayer algorithm applied 
to the two-dimensional Ising and XY models, for various values of the free parameter 
Eq. For Eq — —1 we regain the Metropolis algorithm and for i?o = 1 we regain 
the Wolff algorithm. For — 1 < i5o < 1, we show that the mean size of the clusters 
of (possibly) turned spins initially grows with the linear size of the lattice, L, but 
eventually saturates at a given lattice size L, which depends on Eq. For L > L, the 
Niedermayer algorithm is equivalent to the Metropolis one, i.e, they have the same 
dynamic exponent. For _Bo > 1, the autocorrelation time is always greater than for 
Eq — 1 (Wolff) and, more important, it also grows faster than a power of L. Therefore, 
we show that the best choice of cluster algorithm is the Wolff one, when compared to 
the Nierdermayer generalization. We also obtain the dynamic behavior of the Wolff 
algorithm: although not conclusive, we propose a scaling law for the dependence of the 
autocorrelation time on L. 
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1. Introduction 

Numerical simulations have been widely used in the study of physical systems, specially 
in the last decades. The field of statistical mechanics, among others, has benefited a 
great deal from the use of this technique. In particular, Monte Carlo methods allowed 
for a precise determination of thermodynamic parameters in a variety of models, both 
classical and quantum. Excellent reviews on these methods can be found in Refs. [1] 
and p]. 

In recent years, this field has seen a fast development of new algorithms, which aim 
to make the simulation more efficient, both in time and in memory, as well as broadening 
its application to more complex systems. As examples of these developments, we can 
recall: the calculation of the density of states through fiat histograms, which allows 
obtaining information at any temperature from one single simulation, independent of 
temperature |3]; the use of bitwise operations and storage, which increases by a great 
deal the speed of the update process and saves memory (with the drawback that this 
procedure can be used only with specific models) [1]; and the introduction of cluster 
algorithms, which updates collections of spins, decreasing the autocorrelation time and 
almost eliminating critical slowing down [H El |5l E] . 

In this work we will focus on this last issue. In fact, critical slowing down is a serious 
drawback, which makes simulation of systems at, or near, critical points very inefficient. 
This phenomenon is measured through the scaling of the autocorrelation time, r, with 
the linear size of the lattice, L, assumed to be in the form r ^ L^, for points at the 
critical region. The popular Metropolis algorithm, for example, when applied to the 
Ising model in two dimensions, presents z ~ 2.17 [7J. Algorithms which update clusters 
of spins (the so-called cluster algorithms) have a much lower value of z: this is the case 
for the Swendsen-Wang [5| and Wolff [6] algorithms, for which z is approximately zero 
for the two-dimensional Ising model [8], [9]. 

An alternative (and generalization) to these last two cluster algorithms, the 
Niedermayer algorithm, was introduced some time ago [10] but, to the best of our 
knowledge, has never had his dynamic behavior studied in detail. In this work, we 
calculate the dynamic exponent for this algorithm, applied to the Ising and XY models, 
for some values of the free parameter Eq (see below), in order to determine the best 
choice of this parameter. 

This work is organized as follows: in the next section we present the Niedermayer 
algorithm and relate it to Metropolis' and Wolff's. In Section [3] we review some features 
connected to the autocorrelation time and the dynamic exponent z, in Section H] we 
present and discuss our results, and in the last section we summarize the results. 

2. The Niedermayer algorithm 

The Niedermayer algorithm was introduced some time ago and is an option to Wolff or 
Swendsen-Wang cluster algorithms. The idea is to build clusters of spins and accept 
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their updating as a single entity, hopefully in a more efficient way, when compared to 
these last two algorithms. In this work, we have chosen to build the clusters according to 
the Wolff criterion (they can be constructed according to the Swendsen-Wang rule but 
the results will not differ qualitatively in two dimensions and in higher dimensions Wolff 
algorithm in superior to Swendsen- Wang's). It works as follows, for the Ising model 
(the generalization of this algorithm for the XY model is presented in the Appendix): 
a spin in the lattice is randomly chosen, being the first spin of the cluster. This spin 
is called the seed. First-neighbours of this spin may be considered part of the cluster, 
with a probability 



where K = J/kT, T is the temperature, J is the exchange constant, and Eij is the 
energy between nearest-neighbour spins in unities of J (i.e, Eij = —SiSj; Si,Sj = ±1). 
The free parameter Eq controls the size of the clusters and the acceptance ratio of their 
updating, as seen below. First-neighbours of added spins may be added to the cluster, 
according to the probability given above. Each spin has more than one chance to be 
part of the cluster, since it may have more than one first-neighbour in it. When no 
more spins can be added, all spins in the cluster are fiipped with an acceptance ratio 
A. Assuming that, at the frontier of the cluster there are m bonds linking parallel spins 
and n bonds linking anti-parallel spins, A satisfies: 



where a ^ b represents the possible updating process, from the "old" (a) to the "new" 
(6) state, which differ from the flipping of all spins in the cluster, and 6 — > a represents 
the opposite move. This expression ensures that detailed balance is satisfled [2]. 
Now we must consider three cases: 

(i) for — 1 < i?o < 1, only spins in the same state as the seed may be added to the 
cluster, with probability Padd = 1 — e'^^^^~^^°\ The acceptance ratio (Eq. [2]) cannot 
be chosen to be one always and is given hj A = e-^(i-^o)(m-n)^ if n < m (i.e, 
if the energy increases when the spins in the cluster are flipped), or by A = 1, 
if n > m (i.e, if the energy decreases when the spins in the cluster are flipped). 
If Eq = —1, we obtain the Metropolis algorithm, since only one-spin clusters are 
possible and the acceptance ratio is A = for positive Ai? and 1 otherwise, 
where ^E = 2{m — n) is the difference in energy when the spin is flipped, in units 
of J; 

(ii) for Eq = 1, again only spins in the same state can take part of the cluster, with 
probability Padd = 1 — e~^^. Now, the acceptance ratio can be chosen to be 1, i.e, 
the cluster of parallel spins is always flipped. This is the celebrated Wolff algorithm; 

(iii) for Eq > 1, spins anti-parallel to the seed may be part of the cluster, with probability 
Padd = 1 — e^^^~^°\ while spins in the same state of the seed have a probability 




(1) 




(2) 
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Padd = 1 — e'^^^^^°^ of being added to the cluster. The acceptance ratio is again 
always 1. Note that for Eq ^ 1 nearly all spins will be in the cluster and the 
algorithm will be clearly inefficient (in fact, it will not be ergodic for Eq oo). 
Therefore, we expect that, if the optimal choice of Eq is greater than 1, it will not 
be much greater than this value. 

Our goal here is to do a systematic study of the Niedermayer algorithm, in order 
to establish the optimal value for Eq, at least for the two models addressed in this text. 



3. Autocorrelation time and dynamic exponent 

One possible way to access the dynamic behavior of a numerical algorithm is to measure 

the autocorrelation time, r, of some convenient physical quantity, which is obtained 

from the dependence of the autocorrelation function, p{t), on the time t. Here, time 

is measured in Monte Carlo steps (MCS); one MCS is defined as the attempt to fiip 

spins, where N is the number of spins in the (finite) lattice being simulated (in our 

case, N = L"^, where L is the linear size of the lattice). In fact, a rescaling of the 

time is necessary, when dealing with cluster algorithms [2] and comparing the results 

for different values of Eq. The relation between "time" in MCS, tMcs, and the "time" 

taken to build and possibly flip a cluster, t, is 

< n > . . 

tMCS = t ^ , (3) 

where < n > is the mean size of the clusters. Note that, for Metropolis, < n >= 1 and 
1 MCS is the "time" taken to try to flip spins, as usual. 

In this work, this rescaling has been done and all times are expressed in MCS. The 
function p{t) is deflned as: 

p{t) = J [$(t')- < $ >] [$(t' + t)- < $ >] dt' 

'^t')^t' + t)- >'^]dt', (4) 



where $(t) is some physical quantity. Of course, time is a discrete quantity in the 
simulations; therefore, we have to discretize the previous equation, which leads to [2]: 

pit) = - — - E m'mt'+t)]- 



tmax t 



t'=0 



r, — TT^ E mx E ^(t' + t) (5) 

X'^max ''/ t'=0 t'=0 

The autocorrelation function is expected to behave, as a function of time, as [2] 

Pit) = Ae-'l\ (6) 

at least in its simplest form. It is known that, in some cases, more than one exponential 
term is required [Hj; we will comment on this later. Usually, one can measure r from the 
slope of an adjusted straight line in a semi- log plot of the autocorrelation function versus 
time. However, the autocorrelation function is not well behaved for long times, due to 
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bad statistics (this is evident from Eq. [5l since few "measurements" are available for long 
times). Therefore, one has to choose the region where the straight line will be adjusted 
very carefully and it turns out that the value of r so obtained is strongly dependent 
on this choice. Alternatively, one can integrate p(t), assuming a single exponential 
dependence on (past and forward) time, and obtain: 
1 /■°° p{t) 



T 



-dt, (7) 



2 J-oo p(0) 
with: 

Pit) = e-l*l/^ (8) 
Eq. [71 when discretized, leads to [T2] : 

r--^t^- (9) 

Of course, the sum in Eq. [9] cannot be carried out for large values of t. It has to be 
truncated at some point; we use a cutoff (see Ref.|12] and references therein), defined 
as the value in time where the noise in the data is clearly greater than the signal itself. 
With the value of r obtained as explained above, we made the integral of p{t) / p{0) from 
the value of the cutoff to infinity. A criterion to accept the cutoff is that the value of 
this integral is smaller than the statistical uncertainty in calculating r. Since the value 
we obtain for r is underestimated, this criterion is a safe one. 

Whenever possible, we fitted the autocorrelation time to the expected behavior, 
namely r ~ L^, in the critical region, where z is the dynamic exponent. A point worth 
mentioning is that the autocorrelation function of different quantities may behave in 
different ways. A typical example is shown in Fig. [H where both the magnetization 
and the energy autocorrelation functions are depicted as functions of time, for the 
Niedermayer algorithm with Eq = 0.3 and linear sizes L = 16 (main graph) and L = 256 
(inset). Note the abrupt drop of the magnetization autocorrelation function for small 
times and L = 16. This is an indication that this function is not properly described 
by a single exponential. On the other hand, the energy autocorrelation time follows 
a straight line even for the smallest times. Therefore, we should calculate r from the 
latter, for L = 16, using Eq. [91 However, when L is increased, the picture changes and 
now the magnetization autocorrelation function is well described by a single exponential 
(for small and intermediate values of time), as depicted in the inset of Fig. [H Whenever 
a crossover like this is present, we measure the dynamic exponent from the behavior for 
large values of L and for the function which is well described by a single exponential 
for this range of L, using Eq. [9l But note that, for intermediate values of t, the slopes 
of both curves in Fig. [H (main graph and inset) appear to be the same. However, 
we have already commented on the drawback of calculating r from the slope of the 
autocorrelation function on a semi-log graph. As final notes, we would like to mention 
that we used helical boundary conditions and 20 independent runs (each with a different 
seed for the random number generator) were made for each Eq and L. For each seed, 
at least 4 x 10^ trial flips were made, in order to calculate the autocorrelation functions 



Dynamical heht 




Figure 1. Magnetization and energy autocorrelation functions versus time (in MCS) 
for the Niedermayer algorithm with Eq = 0.3 (see text). The main graph represents 
the behavior for linear size L = 16, while the inset applies to L = 256. 

and their respective autocorrelation times. The values we quote are the average of the 
values obtained for each seed of the random number generator and the uncertainty in r 
is the standard deviation of these 20 values. 

4. Results and Discussion 
4.1. Ising model 

We first present our results for the Ising model and leave to the next subsection the 
discussion of the results for the XY model. 

As already discussed, the case Eq = —1 corresponds to the Metropolis algorithm. 
At the critical temperature, the autocorrelation time scales with L as r ~ L^, with 
z = 2.1665 ± 0.0012 [7]. We have simulated this case only as a test for our algorithm. 
The value we found for z is consistent with the one quoted above and the scaling law 
is obeyed, even for the smallest values of L we simulated. Note also that, for the 
Metropolis algorithm [Eo = —1), it is the magnetization autocorrelation time which is 
well described by a single exponential. 

The first non-trivial value of Eq we simulated was —0.9. In Fig. |2] the 
autocorrelation times for the magnetization is depicted as function of L. We note 
that, for this value of Eq, only the autocorrelation function for the magnetization is 
well described by a single exponential. The initial decay of the corresponding function 
for the energy has an abrupt drop for small times. Therefore, it is not a reliable quantity 
to extract the autocorrelation time from. The value of z was obtained from the curve 
for the magnetization and its value is 2; = 2.16 ± 0.04, which is, within error bars, the 
same value as for the Metropolis algorithm. 

In Fig. [3] the behavior of the mean size of the clusters of spins, < n >, is shown, as 
function of L. For this value of Eq, it seems that < n > does not change with L. We 
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Figure 2. Log-log graphs of magnetization (O ) and energy (□) autocorrelation time 
(in MCS) versus linear size L for the Niedermayer algorithm with Eq = —0.9. The 
quoted value for z is obtained from the slope of an adjusted straight line for the 
magnetization autocorrelation time for L > 16 (see text). The dotted line is just a 
guic' 



<n> 




Figure 3. Mean size of the clusters of possibly flipped spins as function of the linear 
size L for Eq — —0.9. 



will see shortly that in fact it initially grows with L and eventually saturates at some 
value of L, which we call L. 

The overall picture does not change for Eq = —0.5: the magnetization 
autocorrelation function is well described by a single exponential law and the 
autocorrelation time was calculated from it. The dynamic exponent is z = 2.12 ± 0.03, 
still consistent with the Metropolis value (the error bars we quote are all one standard 
deviation; the intersection with the expected value for the Metropolis algorithm, for this 
case, is obtained assuming two standard deviations for the error). Since the picture for 
Eq = —0.5 does not change from the one for Eq = —0.9, we will not depict the graphs 
for the former. 

For Eq = 0, a crossover clearly takes place, as shown in Fig. HI for small L, the 
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Figure 4. Log-log graphs of magnetization (O ) and energy (□) autocorrelation 
time (in MCS) versus linear size L for the Niedermayer algorithm with Eq = 0.0. 
The quoted value for z is obtained from the slope of an adjusted straight line for 
the magnetization autocorrelation time, for values of L beyond the point where the 
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Figure 5. Mean size of the clusters of possibly flipped spins as function of the linear 
size L for Eq = 0.0. 



energy autocorrelation times are larger than their magnetization counterparts, while 
the situation is reversed for larger L (this behavior is more evident for Eq = 0.3; we 
showed the corresponding graph in Fig. [1] above and will comment on it below). The 
value of z is obtained from the slope of an adjusted straight line for the magnetization 
autocorrelation function, for values of L beyond the point where the crossover takes 
place. It reads z = 2.15 ±0.01 in this case, again compatible with the Metropolis value. 
The behavior of < > is shown in Fig. [5l it grows initially with L but eventually 
saturates at L ~ 15. For small values of L it is the autocorrelation function for the 
energy which is well described by a single exponential, while the corresponding function 
for the magnetization shows an abrupt drop for small times. The situation is reversed 
for L > L. 
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This picture is maintained for Eq > 0.0, with the value of L increasing with Eq 
and the crossover taking place at larger and larger values of L. The dynamic exponent 
z is given by 2.16 ± 0.03 and 2.12 =b 0.04 for Eq = 0.3 and 0.5, respectively. Both are 
compatible with the value for the Metropolis algorithm. 

In Fig. [T]we show the change in the behavior of the autocorrelation functions for 
the magnetization and the energy. Note that the crossover mentioned above is connected 
also to the possibility of describing the autocorrelation function by a single exponential: 
this is accomplished by the energy autocorrelation function for small values of L and for 
its magnetization counterpart for larger values of L. 

Finally, for Eq = 0.7 and 0.9 the crossover happens at values of L large enough 
to prevent a reliable estimate of z. It is necessary to go to values of L well above our 
present computational capabilities to be able to extract z from the graphs. 

Nevertheless, the overall trend is well determined: for < Eq < 1, the dynamic 
behavior is the Metropolis' one but this behavior sets in only for large enough L. 
The size of the clusters of turned spins increases with Eq but eventually saturates for 
L = L, where L increases with Eq. For L > L, the relative size of the clusters (i.e, the 
ratio < n > / L"^) decreases and, in this sense, the algorithm is like a single-spin one 
(Metropolis, in our case), explaining the value of its dynamic exponent. Therefore, the 
Wolff algorithm (corresponding to Eq = 1) is still the best choice, when compared to 
the Niedermayer algorithm with ii^o < 1- 

We postpone the discussion of the Wolff algorithm and go to i?o > 1 • In this 
case, spins in different states may be part of the same cluster, although with a smaller 
probability than spins in the same state, and a cluster will always be flipped (see {Hi) 
on page 2). For Eq ^ 1, almost all spins in the finite lattice will take part in the 
cluster and the algorithm will not be optimal (in fact, it won't even be ergodic for 
Eq oo). Therefore, if the Niedermayer algorithm is more efficient than Wolff's, it 
should be for Eq close to 1. We, therefore, studied the cases Eq = 1.1 and 1.05. The 
results are qualitatively equivalent and in Fig. |6]we show both. Note that the growth 
of T with L is faster than a power law for both values of Eq. In the inset, we show the 
corresponding graph for £"0 = 1.1: a crossover is also present but the value where it 
takes place decreases with Eq and for i?o = 1.1 it is not seen. Since the value of the 
autocorrelation time is already greater then for the Wolff algorithm, for a given L, and 
it grows faster than a power law with L, again the optimal algorithm is Wolff's. 

We single out the discussion of the Wolff algorithm {Eq = 1), in order to compare 
with other values of Eq. Although our intention is not to calculate a precise value of 
z for this algorithm, we have adjusted the autocorrelation time for the energy (in this 
case, it is this function which is well described by a single exponential for small values 
of the time) as a function of L for three different functions, namely: 



The first function is the usual scaling law assumed for the autocorrelation function at 




A\nL + C 
A{\nLy + C 



(10) 
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Figure 6. Log-log graphs of magnetization and energy autocorrelation time (in MCS) 
versus linear size L for the Niedermayer algorithm with Eq = 1.05 and = 1.1 (inset). 
In both graphs the autocorrelation time r is plotted as function of the linear size L. 
The dotted line is just a guide to the eye. 



the critical point. We can see in Fig. [7] that there is no indication that the best adjusted 
curve will eventually be a straight line, in a log- log plot. Since previous calculations 
tend to point to a value of z close to zero for the Wolff algorithm in two dimensions, 
one cannot exclude the possibility of a logarithmic dependence, which is the case of the 
second function in the above equation. The fitting is better than for the power law but 
it is not a satisfactory one either. Moreover, it tends to deviate from the data for large 
enough L. The third function is an ad hoc assumption, which proved to be the best fit 
to our data, as can be seen in Fig. [71 The parameters of the function are obtained from 
a non-linear fitting: 

T = AilnLf + C, (11) 

with A = 0.21 ± 0.01, z = 1.50 ± 0.02 e C = 0.47 ± 0.03. We have no theoretical 
explanation for this behavior. The constant C, however, is a finite-size correction. The 
behavior in Eq. [11] is expected to hold true for large enough L and the logarithmic 
dependence makes the scaling region to be reached only for very large values of L. 
For this region, one would expect a simpler law, namely r = A{lnL)^] however, for 
intermediate or small values of L, the constant C acts as a finite-size correction. A 
similar scaling law was found for the exponential relaxation time for the Swendsen- 
Wang algorithm [T3] . 

We also depict the mean size of the clusters of turned spins, < n >, as a function 
of L in Fig. [81 The slope of the straight line is 1.7500 ± 0.0001, which is, as expected 
[2], the value for the ratio 7/z/. Note that, contrarily to what happens for Eq < 1, there 
is no saturation of < n > with L. This seems to explain why Wolff and Niedermayer 
algorithms are in different dynamic universality classes. 
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L 

Figure 7. Log-log graph of the energy autocorrelation time (in MCS) as a function 
of ' — ' r,.. . . . showed. 




L 

Figure 8. Log-log graph of the mean size of turned clusters as function of L. The 
slope is an evaluation of j/v. 

4.2. XY model 

We have applied the Niedermayer algorithm in the study of the dynamic behavior of 
the XY model as well. The generalization of this algorithm to continuous models is 
outlined in the Appendix. Although we have studied three values of Eq, our results are 
conclusive and lead to an overall picture which is analogous to the one for the Ising 
model. We have used the value ksTc/ J = 0.8865 for the transition temperature of the 
two-dimensional XY model. This value is only 0.7% off of the most recent evaluation 
of ksTc/ J for this model [14] . 

The mean size of the clusters of flipped spins for the Wolff algorithm as a function 
of the linear size of the lattice is depicted in Fig. El The slope of the straight line is 
1.7454±0.009, which is slightly different from the expected value for 2 — 77 for this model 
at ksTc/J, 7/4 [15]. In fact, the small discrepancy may be due to the fact that we are 
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Figure 9. Log- log plot of the mean size of the clusters of flipped spins for the XY 
model versus the linear size of the lattice. The slope of the curve just misses the 
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Figure 10. Autocorrelation times for the magnetization (circle) and energy (square) 
for the Wolff algorithm applied to the two-dimensional XY model. 



not using the (unknown) exact value for the transition temperature. 

The autocorrelation times for the magnetization and energy for the Wolff algorithm 
are shown in Fig. [TOl We have not tried to fit the data but it is evident that the 
energy autocorrelation time grows with L slower than a power law. The decrease in the 
magnetization autocorrelation time has been observed previously (in fact, an oscillation 
was observed in an algorithm which mixed Wolff's and Swendsen- Wang's procedures 
but the overall picture is qualitatively similar to ours; see Ref.[T6]). 

We have simulated also the cases Eq = and Eq = —0.5. The mean size of 
clusters of flipped spins saturates and the value of saturation increases with Eq (see 
Fig. [TTl) . Therefore, one expects the same picture as for the Ising model: in particular, 
the dynamic behavior for L large enough is the Metropolis' one. This is confirmed for 
Eq = 0.0 explicitly, where the dynamic exponent measured is z = 1.916 ±0.004 (see Fig. 



Dynamical behav 



13 















p- " 








60 




5 


_l ' 1 ' 1 ' 1 ' 1 ' 


1 ' _ 

o - 






4,8 






<n> 








40 


ja Eg=0.0 


4,6 


- I E=-0.5 








4,4 






20 






& 

' 


1 , " 



I I , I , I , I I I , I , I 

20 40 60 80 100 120 

L 

Figure 11. Mean size of tlie clusters of flipped spins for Eq — —0.5 (inset) and £^o = 




L 



Figure 12. Autocorrelation time for the magnetization and energy as a function of L 
for Eq = —0.5 (main graph) and Eq = (inset) for the two-dimensional XY model. 



[T2l) . Recalling our reasoning for the Ising model for i?o < 1, we can infer that the value 
just quoted for z is an evaluation of the dynamic exponent for the Metropolis algorithm 
applied to the two-dimensional XY model. Since, to the best of our knowledge, there 
is no previous evaluation of z for this model and for the Metropolis algorithm, we have 
made a crude evaluation of z for this case and obtained the value 1.89 ± 0.03, which is 
in agreement, within error bars, with the value we obtained for Eq = 0.0 for large L. 
Clearly, the Wolff algorithm is the most efficient, when compared to the Niedermayer 
algorithm with the two values of Eq quoted above. We have also simulated one example 
of Eq > 1, namely Eq = 1.05. The behavior is qualitatively the same as for the Ising 
model; see Fig. [131 Note that the growth of the autocorrelation time is faster than a 
single power law; in fact, it is well fitted by an exponential. Therefore, also for the XY 
model the best choice is -Eq = 1 (Wolff algorithm), when compared to the Niedermayer 
algorithm. 
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Figure 13. Autocorrelation time for the magnetization and energy as a function of L 
for Eq = 1.05 for the two-dimensional XY model. 



5. Summary 

Wc have studied the dynamic behavior of the Niedermayer algorithm apphed to the 
two-dimensional Ising and XY models. Our main goal is to compare its efficiency with 
the Wolff algorithm's. The latter is a particular case of the Niedermayer algorithm, such 
that a parameter governing the size of the flipped clusters, Eq, assumes the value 1. 

We show that, for — 1 < £'o < 1, the dynamic behavior eventually recovers the 
Metropolis' (Eq = — 1) one. This behavior is linked to the saturation of the mean size 
of the clusters, which happens for all Eq < 1, leading to a decrease of the relative size 
of these clusters when L increases. 

For the Wolff algorithm and the Ising model, we propose an scaling function for 
the autocorrelation time for the magnetization. This choice is an ad hoc one but fits the 
data very well and does not coincide with any function proposed so far in the literature. 
We were not able to make a fitting with the same statistical quahty for the XY model. 

For £'o > 1, the values of the autocorrelation times are greater than those for the 
Wolff algorithm and grow faster than a power law with L. 

Therefore, at least for these two models, the Wolff algorithm is superior to 
Niedermayer 's. 

6. Appendix 

The Hamiltonian for the XY model can be written as: 

n^-j (12) 

<i,j> 

where J is the coupling constant and si is the spin of site i, represented by a unit vector 
in any direction in the xy plane. 
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To start the cluster we randomly choose a preferred direction h and a spin si. This 
spin is the first one of the cluster. Neighbours of si are added to cluster with probability: 



where in the XY model, Eij — —[si • n){s'j ■ n). Note that, ii Eq < 1, only sites that 
have the same component in the direction n as may be added to the cluster. The 
procedure for the construction of the cluster continues as described for the Ising model. 
After the cluster is built the new directions of the spins are given by a reflection with 
respect to axis perpendicular to h. 

The acceptance ratio for XY model is slightly different from the one for the Ising 
model. We cannot define the energy difference (AE) as the number of parallel and 
anti-parallel spins to the cluster. So we must calculate the energy before and after the 
cluster is flipped. In this case we deflne de acceptance ratio, for £'0 < 1, as: 



where a and h have the same meaning as before and /S.E is the difference in energy 
between conflgurations a and 6, in units of J. As we can see, for £'0 = — 1 we regain the 
Metropolis algorithm with A — e~^^^ and for Eq — 1 we regain the Wolff algorithm 
with A = 1 for all clusters. These choices ensure that detailed balance is obeyed. 

The generalization for Eq > 1 is analogous to the one described above and, again, 
spins with different signs for the component along n may also be part of the cluster and 
A = 1 always. 
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